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ABSTRACT 

The  estimation  of  time  delay  and  Doppler  difference  of  a  signal  arriving  at  two 
physically  separated  sensors  is  investigated  in  this  thesis.  Usually,  modified  cross  power 
spectrum  coupled  with  Doppler  compensation  is  used  to  detect  a  common,  passive  signal 
received  at  two  separated  sensors.  Another  successful  approach  uses  the  cross  coherence 
to  achieve  this  goal.  This  thesis  modifies  these  two  techniques  to  model  the  Doppler 
difference  via  an  autoregressive  (AR)  technique.  Analytical  results  are  derived  and  ex- 
perimentally verified  via  a  computer  simulation.  Performance  at  high  and  low  signal  to 
noise  ratios  (  SXRs  )  is  examined. 
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I.      INTRODUCTION 

This  thesis  investigates  the  use  of  autoregressive  (AR)  models  for  estimating  the 
Doppler  difference  and  differential  time  delay  by  processing  of  a  narrow  band  signal 
emitted  from  a  moving  source  and  received  at  two  physically  separated  sensors.  If  the 
signal  is  received  at  two  different  geographical  positions  even  in  the  presence  of  uncor- 
rected noise,  then,  depending  on  the  signal  strength  and  duration,  it  is  possible  to  esti- 
mate the  differential  time  delay. 

Because  the  source  is  moving,  the  signals  that  are  received  at  the  sensors  may  have 
different  frequencies  due  to  the  Doppler  effect.  To  obtain  accurate  differential  time  de- 
lay estimation.  Doppler  difference  compensation  is  usually  required. 

This  compensation  can  be  implemented  by  using  frequency  shifting  of  the  narrow 
band  components  of  the  received  signal.  This  frequency  shift  can  be  obtained  using  a 
Fourier  transform.  In  this  thesis  we  use  an  AR  model  to  detect  the  frequency  shift. 
Using  this  Doppler  compensation,  an  estimate  of  differential  time  delay  can  be  obtained. 
Estimating  the  delay  and  Doppler  using  an  AR  model  can  be  interpreted  as  a  form  of  a 
narrow  band  coherence  procedure,  provided  the  estimations  are  properly  normalized. 

This  thesis  is  arranged  in  six  chapters  and  six  appendices.  Because  the  estimation 
of  the  time  delay  and  Doppler  is  intimately  related  to  the  coherence  between  two  trans- 
formed complex  signals,  an  extensive  investigation  of  coherence  is  given  in  Chapter  II. 
In  Chapter  III.  AR  models,  advantages  of  AR  modeling,  and  AR  model  order  selection 
are  presented.  Chapter  IV  is  devoted  to  the  analysis  and  the  processing  of  noisy  signals 
to  estimate  the  differential  time  delay  and  the  Doppler  difference.  To  estimate  the  dif- 
ferential time  delay,  two  approaches  are  pursued.  AR  model  performance,  Doppler  es- 
timation and  two  types  of  time  delay  estimation  are  examined  in  Chapter  V.  In  the  last 
chapter  some  general  conclusions  of  the  work  carried  out  in  this  thesis  are  presented,  and 
some  suggestions  for  future  investigations  are  given.  Computer  simulation  programs  are 
included  in  Appendices  E  and  F. 


II.      COHERENCE 

A.  DEFINITION 

The  coefficient  of  coherence  between  two  wide  sense  stationary  random  processes 
is  the  normalized  cross  power  spectral  density  function  defined  by  Wiener  [Ref.  1:  p.  12] 
as 

Gxv{f) 

Yxy{/)=       r—-==r  (2-1) 

where  /denotes  the  frequency  [Hz)  , 

Gxy(f)  is  the  cross  power  spectrum  between  x(t)  and  v(/)  , 
GXJJ)  denotes  the  auto  power  spectrum  of  jc(/)  ,  and 
G}Xf)  denotes  the  auto  power  spectrum  of y(i)  . 
Despite  some  confusion  in  the  literature.  Wiener  intended  for  the  coefficient  of  co- 
herence to  be  complex.  This  is  apparent  since  he  discusses  both  the  modulus  and  the 
argument  of  the  coefficient  of  coherence. 

B.  PROPERTY  OF  THE  COHERENCE  FUNCTION 

The  power  spectral  density  matrix  Q(J)  is  positive  semi  definite.  Therefore,  for  two 
random  processes  x  and  y.  we  see  that 


del  [£(/)]  =    del 
For  real  processes  we  have  G  if)  =  G'xv(f)  and  thus 


GJf)     Gxy(j) 
Gyx(J)     Gyy{f)A 


>0  (2.2) 


Gx.(f)Gyy(J)-  lG,,OOl2>0  (2.3) 

where  *  denotes  the  conjugate  of  a  complex  number.   And 

Gxx(J)Gyy(J)>  \Gxy\2  (2.4) 

Furthermore,  GJJ)  and  Gyy{f)  are  nonnegative.  real  functions  of  frequency.  When 
(?„(/)  and  GJf)  are  strictly  positive  definite  (  i.e.,  Gx,(J)GyJf)  >  0  ).  Eq.  (2.4)  can  be  di- 
vided by  Gxx(f)Gyy(J)  without  changing  the  sense  of  the  inequality. 


This  provides  as  an  upper  bound 

k,v(/)l<l         V/  (2.5) 

and  since  the  magnitude  of  any  complex  number  is  greater  than  or  equal  to  zero,  we 
have  the  lower  bound 

0<  \vxy(f)\<\  (2.6) 

The  magnitude  of  the  coherence  function  is  always  between  zero  and  one.  It  is  zero 
if  the  processes  x(t)  and  y(t)  are  uncorrected  and  it  is  equal  to  unity  if  there  exists  a 
linear  relationship  between  x(t)  and y(i)  .  In  order  to  define  the  coherence  it  is  necessary 
that  the  numerator  and  denominator  of  Eq.  (2.1)  are  not  simultaneously  equal  to  zero. 
Coherence  is  not  defined  if  either  auto  spectra  is  zero. 

C.      COHERENCE  ESTIMATION 

If  X ,{/„.)  denotes  the  Fast  Fourier  Transform  (FFT)  of  the  /  th  segment  of  *(«)  at 
frequency/,  ,  then  the  spectral  density  estimates  are  given  by  [Ref.  2:  p.  22] 

v 

/=] 

v 

C,,l/l-)  =  *)  Xtlfk)Y]{fk)  (2.8) 

!=} 
N 

Gyy(fk)  =  ^\yifk)\2  (2-9) 


/=1 


where  a  = 


NTfi  ' 

A  =  number  of  segments, 

T  =  segment  length,  and 
fs  =  sampling  frequency. 

Finallv  the  coherence  estimation  is 


to-  ^^ 


\    Gxx(f^Gyy(fk) 
,V 


/E|a™|2E|}^)|2 

D.      COHERENCE  OF  NARROW  BAND  SIGNALS  WITH  DIFFERENTIAL  TIME 
DELAY  AND  DIFFERENTIAL  DOPPLER 

The  output  of  the  band  pass  filters  (BPF1  and  BPF2)  in  Figure  1  are  denoted  by 
X/(fk)  and  Y,{J])  respectively.  Each  term  represents  the  Fourier  transform  of  the  corre- 
sponding time  series  evaluated  at  frequency/,  and  time  /.  For  narrow  band  signals  a 
Doppler  shift  corresponds  to  a  frequency  shift.  If  a  signal  arrives  at  the  two  sensors 
having  a  differential  Doppler  shift  as  well  as  a  differential  time  delay,  then  we  see  that 
a  frequency  compensation  and  time  delay  compensation  by  the  appropriate  values  tend 
to  line  up  the  signals  in  frequency  and  time.  This  is  accomplished  by  using  an  additional 
Fourier  transform  in  channel- 1  of  the  processor  and  a  time  delay  in  channel-2  of  the 
processor.    Mathematically,  this  can  be  expressed  as 


>■(/"*)=      p^—        v  (2-11) 

V    /=]  /=1 

Comparing  this  with  Eq.  (2.10),  we  see  that  we  generalized  the  coherence  concept.  We 
also  note  that  the  implementation  resembles  a  correlator,  where  one  of  the  signals  is 
frequency  compensated  and  the  time  delay  corresponds  to  the  delay  operation  of  the 

correlator. 
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Figure   1.        Coherence  estimation  block  diagram. 


If  the  Figure  1    is  redrawn  as  in  Figure  2,  then  it  can  be  interpreted  as  an  FFT 
implementation  as  shown  in  Figure  3. 
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Figure  2.       Coherence  estimation  block  diagram  (reinterpretation). 
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Figure  3.       Coherence  estimation  block  diagram  using  the  FFT. 

But  this  FFT  has  a  poor  resolution.    To  obtain  a  better  resolution,  an  AR  model 
is  desirable.    This  approach  is  shown  in  Figure  4. 
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Figure  4.       Coherence  estimation  block  diagram  using  an  AR  model. 


Figure  1  and  Figure  4  represent  two  different  schemes  to  compute  the  coherence 
function.  Differential  time  delay  and  Doppler  difference  can  be  estimated  by  using  AR 
models  in  the  modeling  of  the  coherence  function. 


III.      AUTOREGRESSIVE  (AR)  MODELS 

A.      AR  MODELING 

If  the  following  difference  equation  is  satisfied,  the  resulting  structure  is  called  an 
AR  model  of  order  p.   [Ref.  3:  p.  177] 


in)  =  -  /  PiAn  -k)  +  u{n) 


k=\ 


(3.1) 


where  x{n)  =  the  signal  at  instant  n, 
u{n)  =  the  white  noise  driver,  and 
ak  =  the  k  th  coefllcient  of  the  AR  model. 

A  realization  of  the  AR  model  is  illustrated  in  Figure  5. 
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Figure  5.       AR  filter  of  order  p. 


B.  ADVANTAGES  OF  THE  AR  MODEL 

The  motivation  for  parametric  modeling  of  random  processes  is  the  ability  to  obtain 
better  spectral  estimates  based  upon  the  model  than  estimates  produced  by  classical 
spectral  estimation.  Both  the  periodogram  and  correlation  methods  can  be  used  to  yield 
Power  Spectral  Density  (PSD)  estimates.  The  unavailable  data  or  autocorrelation  se- 
quence (ACS)  values  outside  a  given  window  are  assumed  to  be  zero.  This  kind  of  an 
unrealistic  assumption  leads  to  distortions  in  the  spectral  estimate. 

The  advantages  of  the  AR  approach  are 

1.  AR  spectra  tend  to  have  sharp  peaks,  a  desirable  feature  of  high  resolution  spectral 

estimators. 

2.  AR  parameter  can  be  obtained  as  solutions  to  linear  equations. 

C.  POWER  SPECTRAL  DENSITY 

Eq.  (3.1)  can  be  rewritten  as  follows 

p 
x(n)  =  -  /  aIcx{n  -  k)  +  u{n) 


=  Thku(n-k)  (3.2) 


k=0 


where  hk  =  the  causal  filter  impulse  response. 
Let  us  take  the  Z  -transform  of  Eq(3.2). 


X{z)  =  H{z)U(z)  (3.3) 


Rewriting  Eq.  (3.1)  as 


Yf,r<n-k)  =  u{n)  (3.4) 

fe=0 

where  a0  =  1 ,  and  taking  the  Z  -transform  gives 

A(z)X(z)  =  U{z)  (3.5) 

Eq.  (3.3)  and  Eq.  (3.5)  can  be  solved  to  obtain  H(z) 


m-jfi  (3.6) 

and  hence 

-V(z)  =  -^  (3.7) 

The  Z  -  transform  of  the  output  sequence  {x(n)}  is  related  to  the  Z  -  transformation 
of  the  input  random  process  u{ri)  by  [Ref.  4:  p.  56] 

U\z)U{±-) 

pxx(z)  =  x{z)x\  4- ) = Hi— = pJ&  — lt_         (3-8) 

A(z)A(-±r)  A(z)A(-±r) 


The  AR  power  spectral  density  is  obtained  by  substituting  z  =  ennfT  into  Eq.  (3.8) 
and  scaling  it  by  the  interval  T. 

Par(J)  =  Tpm  j-^--  =  TPw  '  (3.9) 

Ml/)  I  ep  (f)aa   ep(f) 

where  er  =  [  l.e-*"*7 e-M»T3H 

a  =  f  1.  a,  c7, a,,2r 

pa  =  variance  of  driving  sequence. 

D.     BURG'S  ALGORITHM 

In  practice,  the  autocorrelation  is  usually  not  available,  so  one  must  make  an  AR 
spectral  estimation  based  on  the  available  data  samples.  The  simplest  procedure  to  ob- 
tain an  AR  spectral  estimate  from  data  samples  would  be  to  produce  estimates  of  the 
autocorrelation  sequence  from  the  data.  These  autocorrelation  estimates  would  be  used 
in  lieu  of  the  true  autocorrelation  sequence  in  the  YULE  -  WALKER  equations  to  yield 
the  AR  coefficients.  However  better  results  are  obtained,  particular}'  for  short  data 
segments,  by  algorithms  that  obtain  the  AR  model  parameters  directly  from  the  data, 
without  explicitly  estimating  the  autocorrelation  function. 

The  Levinson  recursive  solution  to  the  YULE  -  WALKER  equations  relates  the 
/>th  order  AR  parameters  to  the  p  —  1  th  order  parameters  as  given  by  [Ref.  3:  p.  21 1  ] 

Op{n)  =  ap_x{n)  +  kpap_x{p  -  n)  (3.10) 


For  n  =  1  to  n=p—l,  the  reflection  coefficients  kp  can  be  found  by  using  the 
known  autocorrelation  function  for  lag  0  to  p  —  1.  So  we  have 


P-\ 


)j*p-\(n>jJP  ~  n) 


kp  =  ap(P)=       "=Q  (3.11) 

The  recursion  for  the  driving  white  noise  variance  is  given  by 

PP  =  Pp-,(\-  \kp\2)  (3.12) 

where  p0  =  r„(0)  . 

But  the  ACS  is  not  available,  hence  we  can  not  calculate  the  reflection  coefficients. 
The  Burg  algorithm  provides  an  estimate  of  the  reflection  coefficients  which  in  turm  are 
obtained  through  a  least  squares  criterion. 

The  forward  linear  prediction  error  is  given  by 


/jn)  =  x(  n)+    }  (/p{m)x{n  -  m)  (3.1 3 ) 

while  the  backward  linear  prediction  error  is  given  by 

p 
4^  =  -v("  ~  P)  +  /   4>*{™)x{n  +  m  ~  P)  0- 14> 

Substitution  of  Eq.  (3.10)  into  Eq.  (3.13)  and  Eq.  (3.14)  yields  the  recursive  relationships 

/Jin)  =  efp_,{n)  +  k/p_x{n  -  1)  (3.15) 

ebp{n)  =  ebp_x{n-\)  +  k/p_x(n)  (3.16) 

At  each  order  p,  the  arithmetic  mean  of  the  forward  and  backward  linear  prediction 
error  power  (  sample  prediction  error  variance)  is  given  by 

10 


fb    1 
Pp=  — 


n=p+]  '      n=p+] 


(3.17) 


This  expression  is  minimized,  subject  to  the  recursion  given  by  Eq.  (3.15)  and  Eq. 
(3.16).  Thus,  p-f  is  a  function  of  single  parameter,  namely  the  complex  valued  reflection 
coefficient  kp.   Setting  the  complex  derivative  of  Eq.  (3.17)  to  zero 


>rf 


6Re(U      J  clm(kp) 


-,   fb 

0  (3. IS) 


and  solving  for  k  yields 


-2    )    fi£_i(/0£i(« 


*„  =  — ^ — (3.19) 

»=p+i  «=/-'+ i 

The  estimation  of  the  reflection  coefficient  represents  the  HARMONIC  mean  be- 
tween the  forward  and  backward  partial  correlation  coefficient,  where 
e{{ri)  =  eftji)  =  x(n)  and  N    =  number  of  data  points. 

E.      FINAL  PREDICTION  ERROR  (FPE)  CRITERION 

Because  the  best  choice  of  filter  order  is  not  generally  known  a  priori,  it  is  usually 
necessary  in  practice  to  postulate  several  model  orders.  FPE  is  a  kind  of  criterion  which 
was  provided  by  Akaike  [Ref.  3:  p.  230J.  This  criterion  selects  the  order  of  the  AR 
process  so  that  the  average  error  variance  equals  the  sum  of  the  power  in  the  unpre- 
dictable part  of  the  process  and  of  a  quantity  representing  the  inaccuracies  in  estimating 
the  AR  parameters.   The  FPE  for  an  AR  process  is  defined  as 

pnn  a      A     A' +(/>+!) 

FPE(p)  =  pp  x_{p+l)  (,.20) 
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where  A"  is  the  number  of  data  samples, 

p  is  order  of  the  filter,  and 

pp  is  the  estimated  white  noise  variance  when  using  a  p  th  order  filter. 
The  order  p  selected  is  the  one  for  which  the  FPE  is  minimum. 
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IV.      DOPPLER  AND  DIFFERENTIAL  TIME  DELAY  ESTIMATION 

A.      DOPPLER  ESTIMATION 

Let  x{t)  and  y(t)  be  the  signals  received  by  two  sensors 

x(t)  =  cos{2n(f+  a,)*}  (4.1) 

y(t)  =  cos{2n(f+a2)t}  (4.2) 

where /is  the  carrier  frequency. 

a{  ,  c2  are  Doppler  shifts,  and 
i  is  the  time  variable. 
Let  fs  be  the  sampling  frequency  which  satisfies  the  Xyquist  theorem.  For  con- 
venience, in  all  derivations  and  simulations  a  sampling  frequency  of  64  Hz  is  used,  to- 
gether with  band  pass  filter  width  of  1  Hz.  We  assume  that  |a,|  <  0.5  (  /=1,2  )  . 
which  implies  that  the  signals  stay  in  the  band  pass  filter  regions  of  their  respective  band 
pass  filters  regardless  of  any  Doppler  shift.  Note,  these  values  can  be  modified  to  arbi- 
trary sampling  rates  and  pass  band  regions.    The  sampling  rate  is  given  by 

fs>2f+l>2(f+v.i)  (    i=lt2    )  (4.3) 

The  phase  at  instant  k  and  sampling  time  interval  T  are  given  by 


(4.4) 

(4.5) 
(4.6) 


x^k 

2n(f+  a 

{)k 

fs 

y<f>k  : 

2n(f+  a. 

2)k 

fs 

H 

Let  x(ri)  denote  the  sampled  analog  signal  x(t)  |  ^  then 


n 


x(n  +  m)  =  cos{2n{f+  aj-f  +  x<f>m}  (4.7) 

Js 


13 


y{n  +  m)  =  cos{2n(f+  a2)  jr  +  y4>m)  (4.8) 

Js 

The  derivations  of  Eq.  (4.7)  and  Eq.  (4.8)  are  given  in  Appendix  A. 

These  signals.  {x(n)}  and  [y{ri)},  are  processed  to  detect  the  Doppler  difference  . 
Let  BPF1  be  the  band  pass  filter  centered  at /J  and  BPF2  be  the  band  pass  filter  centered 
at/,  (vvhere/;  might  equal/  ).  N  data  points  are  processed  in  the  band  pass  filters,  to 
generate  one  output.  The  inputs  to  BPF1  and  BPF2  at  time  m  are  vectors  Xm{ri)  and 
Yjn)  ?  respectively.  The  size  of  the  input  vector  to  the  filter  is  the  number  of  data  points 
taken  during  1  second  (i.e..  the  number  of  points  processed  in  the  filter  =  A"=/  ).  The 
input  vectors  are  denoted  by 

XJn)  =  [-v(m).  x(m  +  1) x(m  +  N  -  1)]  (4.9) 

Yjn)  =  [y(m),y(m  +  1). ,y(m  +  X-  1)]  (4.10) 

The  filtering  is  performed  using  FFTs  ,  where  successive  FFT  outputs  are  generated  at 
the  input  data  rate.    The  BPF2  output  is  conjugated. 

To  avoid  the  complexity  the  following  four  complex  constant  variables  are  defined. 

:V-1 

Ax=  ycos{2^/+a1)7-}^"y2"/^r  (4.11) 


.V-l 

Bx  =  Tsm{2n(f+  a,)  f  }e~j2r:/f  (4. 12) 


,V-1 

Ay  =  Vcos{2tt(/+  aj  f  }e/2lT/f  (4.13) 


,V-l 


By=  ysin{2rr(/+a2)^-}e/2r:/f  (4. 14) 

At  instant  m  the  complex  output  signal  of  BPF1  can  be  calculated  as  follows 
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XJf  =  V cos{2tt(/'+  ai)  7-  +  ^J*"**'  -v 

■V- 

=  /[  cos{2re(f +  aj)y- }  cos  >„,  -  sin{2^(/+  aj  -~  }  sin  ^J^"72^ 

A'-l  A-'-l 

Y"^  ^;         _ji-JL.  V^  f?         _;->_JL 

=  cosx</)w2^co^27r^+ai)7r^    " N  -smx<Pm2^s'm^27:(f+a^Y^e     -v 

(4.15) 

=  ^  cos  ^  -  £x  sin  X0W 

^«*«  4.  e-^<Pm  e/\Om  _  e-jx4>m 

=  Ax  ^  ~  Bx  ~, 

4/ 

Ax  +JBx    J,on.    ,     Ax  ~JBx      -jx<pm 
2  2 

^.r     J^x    j2<p m    .     "-x     J^x      —flq> r-1— m 

= 5 e   '      /,  + e      ■      J: 

The  complex  signal  XJf  contains  two  frequencies  with  different  amplitudes. 

To  understand  the  character  of  XJf)  ,  it  is  important  to  evaluate  which  is  the 
dominant  term  in  the  above  expression. 


.v-i 
Ax  +JBX  =  )[cos{2n(f+  a])f}+j  sm{2n(f+  aj  f  }]e~^l 

i — l  Js  Js 

A'-l 

n        ,.    ,  ?! 


^+*ne**fN 


«_=o  (4.16) 

1  -  ^ 

l-^'l 
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A'-l 


Ax  -jBx  =  2JL  cos{2n(f+  Kl)  jr  }  -J  sm{2r:(/  +  a,)  j  }>_y'2r:/~ 

n=0 

A-l 


n        ■-,   f  n 


-Jl^+afiJL.      nifJL. 


«=o  (4.17) 

A-l 

1  -  c-j1^ 

From  Appendix  C 

l^+Al  >  1 4c -A  I  (4.18) 

Therefore #*e«    is  the  dominant  term  of  XJJ). 

In  similiar  way.  the  complex  output  signal  of  BPF2  can  be  calculated  as  follows 


v-i 


Ym[f)  =     "cos{27r(/+  a2)  -*-  +  y4>m}^f  N 


n=0 

;V-1 


=  ^[  cos{2K(f+  «2)  J  }  cos^  -  sm{2n(f+  a2)  -*-  }  sin^]^  ,v 

/;=0 

A'-l  A-l 

cos{2tc(/"+  a2)  —  K  "  a-  -  sinv(/)w2^sin{27rC^-fa2)7-}^      '  v 

S  (4.19) 

=  .-i,  cos  v(f)m  -  By  sin  v<£, 


";' 


>  ->  ">  2/ 


A  +JBV 


I     .  Av  -jBv       .  . 


Av+jBv     ,.,  y+»z)  A,-jBv       ,. ,  if+  »2> 

=     ;     e/2<3-r~"+  -1 ;    g-^^rw 

The  complex  signal  }'^(/)  contains  two  frequencies  with  different  amplitudes. 
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Even  though  Y'Jf)  is  similiar  in  form  to  XJJ),  it  is  not  obvious  which  term  is  dom- 
inant. 


Ay  +JBy  =  2^[  cos{2;r(f +  a2)  j  }  +j  sm{2n(f+  a2)  j-  Ue12^ 


n=0 

A-l 
\ — > 

t 

^ 

•(f+^J^ 

rr/- 

A' 

L - 

/7=0 

s-\ 

I 

e* 

:(2/+a,)-^- 

«=0 

1 

1 

—  C 

/2r(2/+  *2) 

1  - 

V2 

^(2/+a2)Jr 

1  - 

/2-r^ 

1  - 

V2 

r(2/+  a2)  -L 

4,  -jBy  =  V [  cos{2-(/  +  a2)  -f  }  -j  sin{2^+  a2)  -f  }]^'2rr/  .v 


«=0 
A  - 1 


n=0 
A-l 


^"^^2-y 


«=0 

1  -e 


l-e-J^2- 


(4.20) 


(4.21) 


From  Appendix  C 

MV.-7A.|  >  |  A,  +JBy  |  (4.22) 

Therefore     — — e^*™    is  the  dominant  term  of  Y'JJ). 
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Since  the  two  output  sequences  {XJJ)}  and  { Y'Jf)}  have  two  sinusoidal  components 
each,  their  product  {XJj)TJJ)}  has  four  sinusoidal  components.  The  four  frequencies 
are 

{y.,-a2)  (a2-«P  (2/+ «,  +  «2)  (2/+ftl+«2) 

(ox  =  2n ,  co2  =  2- ,  o>3  =  In ,  a>4  =  -2rr 

Js  Js  Js  Js 

with  <u,  the  sinusoid  with  the  largest  amplitude  and  co2  the  sinusoid  with  the  smallest 
amplitude. 

The  product  of  the  output  of  the  filters  is  given  by 

Ax+jBx     ■  ■     Ay—jBv    e-jy4>m 
XJJ)  Ym(f)  =     x  2        Jx°m  — — 2 —  +f(ot'-  **  *^™'  y^rn) 

(Ax+jBx){Ay-jBy)     ...        .  ,  (4'23) 

=  — "7" "  e^--^  +F(al8  a2,  ^  ,0W) 

where  F(al5  a;.  „$ms  v(/>  J  represents  the  three  low  amplitude  frequency  terms. 

When  using  the  AR  model  as  described  in  chapter  III.  at  a  high  SXR  (  i.e.. 
SXR  ->  00  )  .  a  4th  order  AR  model  detects  all  of  the  four  frequencies  given  above. 
When  the  SXR  is  low  (i.e..  20  dB  ).  only  one  dominant  sinusoid  is  detected  with  fre- 
quency 

Intf+OLi)        2n{f+<x2) 


=  2n 


fs 

Since  we  can  detect  a>,  and  know  the  value  of/  ,  a,  —  a2  can  be  estimated  by  using  an 
AR  model. 

B.      DIFFERENTIAL  TIME  DELAY  ESTIMATION 

To  detect  and  localize  a  target,  it  is  important  to  find  the  differential  time  delay  of 
signals  arriving  at  two  sensors.  Let  us  assume  that  signals  at  the  two  sensors  are  as 
given  in  Figure  6  (  i.e..  zero  differential  time  delay  )  For  the  remainder  of  the  figures,  the 
time  axis  is  scaled  to  be  64  points  per  second. 


% 

X 

1           1  1 1  j  1     ] 

I   I    I   I       III 

1     1 

* 

o 

X 

S*5 
o 

- 

:       111  || 

fill  lis 

III 

z 
1^0 

I  fni 

ll 

X 

III     HI   III   |[| 

ill!  lilt  il  til  iiii il 

1  I 

in 

■a 

I  III! 

X 

1             11 

i  1 1  i  [  1 1  {  i 

f  1 

T 

1                                  1                                  1     1 

i ' 

i                       i                       i 

-100                                           0 

100 

200 

TIME 

% 

- 

1    1 

Il         1   1     ll 

X 

til  il  i! 

111      Mill! 

^" 

1  jll  III  Ij 

Ill   II     1     11     1 '  i 

o 

III    In      it    111    IL  1 1 

| 

Il  I1 1 ;!!!  1 

\ 

O 

'■  mm  pp 

IliillPsIt  <!i|i;      ' 
t.Ub       1  i|!|liSi;ili.tf: 

ill      I1  ! 
jj  ll  1  I  1  jll 

1 

to  o 

z 

X 

;!!l 

m 

■a 

!i  I'M 

1    l! 

X 

1 1 1 1 

1                        f 

7 

1               i            , ,  i 

1 

!                                         |                                         1 

-100                                           0 

100 

2C0 

TIME 

Figure  6.        Receiving  signals  at  sensor  1,2  (SNR=  100,  no  differential  delay). 

We  can  estimate  the  differential  time  delay  and  differential  Doppler  using  two  dif- 
ferent approaches.  The  first  approach  uses  the  cross  power  spectrum  while  the  second 
approach  uses  the  coherence. 

1.  Differential  time  delay  and  differential  Doppler  estimation  using  the  cross  power 
spectrum. 

Figure  4  shows  the  block  diagram  of  the  coherence  estimator  using  AR  mod- 
eling. In  this  figure,  the  output  of  the  FFT  can  be  interpreted  as  the  cross  power  spec- 
trum. Using  this  cross  power  spectrum,  we  can  estimate  the  differential  time  delay  and 
the  differential  Doppler.  But  when  we  use  the  highest  peak  of  the  cross  power  spectrum, 
the  peak  is  somtimes  not  detected  at  the  proper  time  delay  nor  at  the  proper  Doppler 
shift.  We  will  show  a  special  case  in  which  the  peak  of  the  power  spectrum  is  not  de- 
tected at  the  proper  time  delay  nor  at  the  proper  Doppler  values. 


19 


a.      Special  case. 

For  our  test  case  the  data  duration  is  six  seconds,  two  seconds  of  the  noisy 
signal  and  four  seconds  of  noise  only.  Linear  transformation  of  this  data  leads  to  one 
of  three  types  of  outputs.  The  first  type  represents  full  information,  while  the  second 
type  represents  partial  information.  The  third  type  represents  a  noise  only  condition. 
Generally  when  two  signals  are  lined  up  in  time,  the  AR  model  should  give  the  highest 
output  power.  But  this  is  wrong  in  some  cases.  In  the  full  information  case,  the  mag- 
nitude of  transformed  signal  is  high  and  constant.  For  some  reduced  information  case, 
even  though  the  magnitude  of  the  transformed  signal  decreases,  it  still  may  have  large 
magnitude  of  spectral  components.   This  phenomenon  can  be  explained  as  follows. 

Define  two  functions  depending  on  k 


fV-l-A 


k 


Av=    )     cos(2rr/;-7-)e/'2r/A-  (4.25) 


fco  *  & 


.V-1-/C 

D     _ 


Bv=    Y  sin^-^y2*7!!'  (4.26) 

,7^0  J* 


where      f2  —f+  &2. 

k  denotes  the  number  of  lost  data  at  BPF2  input  vector. 

We  assume  that  we  know  when  the  BPF1  signal  starts.  If  this  information 
is  not  available,  we  need  to  examine  the  signal  at  the  output  of  BPF1.  to  obtain  a  can- 
didate time  frame.   Consider  the  input  of  AR  model  as  shown  in  Figure  7. 

The  input  data  size  to  the  AR  model  is  the  number  of  linearly  transformed 
data  points  during  a  given  period  (i.e.,  input  size  is  N=fs  )•  The  BPF1  output  represents 
full  information  (i.e..  each  output  element  is  produced  from  the  signal  information 
points).   Therefore,  as  shown  in  the  previous  section,  the  dominant  output  sequence  is 

kJ\    kJ**\    kJ^\    kJ^\ ,    kJ*"~x  (4.27) 


where 


1       1  _  Jlm\ 
*i-y  t  (4-28) 

-    i  _  (J1^!  x 
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Figure  7.       AR  model  and  its  driver  source. 


(I)  0  second  delay  case.  We  already  defined  the  input  vector  XJn)  and 
Ym{")  to  BPF1  and  BPF2.  This  vector  can  be  interpreted  as  a  time  snap  shot.  Each 
output  of  the  BPF  requires  N  input  points.  If  we  require  N  output  points  from  the  BPF, 
even  using  maximum  overlap  in  the  processing,  we  require  at  least  2A7  —  1  points  at  the 
BPF  input.  Figure  8  shows  the  2Ar  input  to  BPF1  and  BPF2.  Both  2  A"  input  points 
contain  the  signal  and  some  noise.  As  shown  in  Figure  8,  in  the  zero  delay  case,  the  two 
input  vectors  Xm{n)  and  YJn)  (  0  <  m  <  A*  )  provide  full  information.  So  the  lin- 
early transformed  output  Xm(fk)  and  Y'n(fk)  also  represent  full  information. 
The  dominant  part  of  the  BPF2  output  is  the  sequence 


,-/° 


-Mi 


V    >  *2«  '  .  V  '  .  M 


-Jy<t>3 


~U4>i 


.,    k2e  J 


iV-I 


K29) 


where 


*2=, 


1       \-e 


-J2m3 


l-e-***!? 


(4.30) 
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Figure  8.       Tuo  signals.  0  second  delay  at  the  sensors. 

The  dominant  input  sequence  to  the  AR  model  is  given  by 

{klk2e/{'4,k~^k)    .k  =  0.1. 2 -V-1}  (4.31) 

where  the  magnitude  is  kxk2 

^2 j  -I  second  delay  case.  Figure  9  shows  the  2.V  input  to  BPF1  and 
BPF2.  The  2\  input  data  points  of  BPF1  contain  the  signal  and  noise.  During  the  first 
.V  points  the  input  to  BPF2  contains  the  signal  plus  noise  while  during  the  second  A* 
points  only  noise  is  present.  The  input  vector  XJn)  of  BPF1  contains  full  information, 
but  input  vector  Ym{n)  of  BPF2  does  not  contain  full  information.  In  the  Y0{n)  case, 
every  element  of  Y0{n)  is  full  information.  But  when  k  is  not  equal  to  zero  then  Yk(n) 
consists  of  .V—  k  information  elements  and  k  noise  elements. 
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Figure  9. 


Two  signals.  -1  second  delay  at  the  sensors. 


The  output  from  BPF2  represents  partial  information  while  the  out- 
put from  BPF1  represents  full  information.   The  BPF2  output  sequence  is  given  by    . 


{ 


/A  +//A 


&*  + 


kAv-jkBv     _ 


e  -Jy<Pk 


.     A-  =  0.1 V-  1    } 


(4.32) 


The  derivation  of  Eq.  (4.32)  is  given  in  Appendix  B. 

Assuming  that  — — ^ — -  e~h^  is  the  dominant  term,  which  is  a  valid  assumption  for  our 
test  case  with  fs  =  MHz,  f=  21Hz,  a,  =  0.4Hz  and  a2  =  0.001Hz,  then  the  input  se- 
quence to  the  AR  model  is  dominated  by 


{*,  kAy   JkBy  ^-^J      A  =  0,1. 2 A'-l} 


(4.33) 


with  A',  is  defined  earlier  (  Eq.  (4.2S)  ). 
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The  magnitude  of  dominant  input  sequence  is  obtained  by  examining 

\-\-k 

kAy  -JkBy  =    V  [  co${2r,(f+  <x2)}  -j  sin{2n(f+  7.2)}]^'2^f 

n=0 

A-l-A- 


n     .,    ,  n 


n=(i 


K-\-k 


(4.34) 


e-P^2- 


n=0 

„     X-k 
l-c-^  —  'i 


l-e-J^2~ 


When  we  compare  the  magnitude  of  the  sequence  in  Eq.  (4.27)  and  the  magnitude  of  the 
dominant  term  in  Eq.  (4.32),  the  magnitude  of  I  A,  I  >  \ kAy —jkBy\  j2  where 

,         1       1  -  e~J2^- 
K-) 


2  -r-y  -L 

1  -  e  j2"*2  .v 


and 

h'-k 


k 


■\  -JkA         l     1  -  e~nr-  '  v   a2 


2  l-e-Jl«*2JT 


JO  t  -/2^33  k 


(4.35) 


">  »        l  7  ■■>        i 

1  —  e         2  A  \  —  e         2  A 

a2 

We  note  that  Eq.  (4.35)  has  two  frequencies  (  i.e.,  0  and  -rr  2tc  ). 


The  magnitudes  of  both  terms  are 


1  -  r'2"2y 


Now,  note  that  as  Ions  as  I  a,  I  <  — ,  then  the  magnitudes  of  both  terms  in  Eq.  (4.35) 

6 

are  lareer  than    I  k-,  I 


Due  to  the  symmetry  a  positive  delay  results  in  a  similar  reponse  to  a  negative  delay,  and 
hence  Eq.  (4.35)  holds  also  for  delay  of  +  1  second.  Equation  (4.35)  can  be  interpreted 
as  having  large  responses  at  shifted  Doppler  values  at  a  delay  of  1  second  and  a  delay 
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of  —  1  second.  Proper  delay  can  not  be  established  nor  can  the  Doppler  value  be  estab- 
lished. This  cross  power  spectrum  algorithm  will  not  give  good  information  about  either 
Doppler  difference  or  time  delay. 

b.      Modified  cross  power  spectrum. 

To  detect  the  signal,  input  signals  to  the  AR  model  should  be  preprocessed 
as  shown  Figure  10.  If  the  BPF2  output  signal  is  magnitude  normalized,  then  this  nor- 
malized signal  retains  good  phase  information.  When  we  normalize  the  BPF2  output, 
one  of  two  conditions  can  occur.  The  first  condition  (  i.e.,  synchronized  )  leads  to  the 
magnitude  normalization  of  Eq.  (4.29)  and  provides  accurate  frequency  and  delay  esti- 
mation. The  second  condition  occurrs  when  the  signal  in  channel-2  is  not  lined  up  in 
time  with  the  signal  in  channel- 1  (  i.e.,  during  the  delay  search  ).  Under  this  condition 
smaller  peaks  in  the  time  Doppler  plane  appear  at  incorrect  values  of  time  delay  and 
Doppler.  Above  70  dB  SXR,  the  correct  peak  is  the  dominant  one  and  allows  proper 
estimation.  Information  from  contour  plots  can  be  used  at  values  of  SNR  between  70 
and  10  dB  (  see  Appendix  D  ). 
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x,(fk) 
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AR 
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BPF2 
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y(n)      ' 

1 OM 

■< 

Y,'JU 

Figure   10.       Modified  cross  power  spectrum  block  diagram. 

2.      Differential  time  delay  and  differential  Doppler  estimation  using  the  coherence. 

Another  way  detecting  the  signal  is  normalization  of  the  AR  output  with  the 
squared  root  of  the  product  of  PXJJ)  (  auto  spectrum  of*  )  and  Pyy(f)  (auto  spectrum  of 
y  ).   This  can  be  done  using  two  additional  AR  models  as  shown  in  Figure  1 1. 
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Figure   11.       AR  model  of  coherence. 

Using  the  AR  approach,     pxx  ,  pxyxy  ,  pyv  ,  and  AR  coefficients  can  be  ob- 
tained.   Therefore  we  get 


p^in  = ,  Tpx:,2 


<lxxV) 


(4.36) 


Tpyy 

Pyyif)  = ^T 


(^.37) 


'  xyxvv) 


TP 


xy  xy 


^xyxvi/) 


(4.38) 
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where  P„(/)  is  the  power  spectrum  of  A', , 
Pyy{J)  is  the  power  spectrum  of  T,_d  , 
PXVXJJ)  is  the  power  spectrum  of  X,T^d  , 
A„(f)  is  the  AR  power  spectral  density  of  A")  , 
Ayy{f)  is  the  AR  power  spectral  density  of  Y]_d . 
dxyxy{f)  is  ^e  AR  power  spectral  density  of  X,T,^  cross  term. 
pxx  is  the  driving  noise  variance  of  A'-  . 
pyy  is  the  driving  noise  variance  of  Y'^  , 
pxyxy  is  the  driving  noise  variance  of  Xtf]^  .  and 

r-Jr. 

f 

We  obtain  a  coherence  estimate  from  the  three  power  spectra,  by  assuming  that 
|/U/)la*lAwWI    • 


(4.39) 


A2      1  ^(/)  1 2 

1  Pxyxyif)  1               1 

pw     \Axx{j)Ayy{f)\1 

~     l^xx(/)^(/)l           T 

PxxPyy           \  Axyxy{f)  | 

1   1       ^W 

\Axx{f)Avv(J)\2 

T    PxxPyy  \Axyxy{f)\ 


(4.40) 


This  approach  provides  good  information  about  the  differential  time  delay  and 
differential  Doppler  down  to  SXRs  of  20  dB. 
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V.      RESULTS 

This  chapter  presents  graphical  results  obtained  by  applying  the  analytical  results 
of  the  previous  chapters  to  specific  examples  and  carrying  out  the  required  computations 

on  the  computer. 

A.      AR  .MODEL 

This  section  provides  the  AR  model  performance  tests.  As  we  discussed  in  chapter 
III,  Burg's  algorithm  is  applied  to  this  AR  model.  The  FPE  criterion  is  use  to  find  the 
AR  model  order.  Figure  12  is  the  power  spectrum  of  two  test  signals  s'mico^Tn)  and 
cos(fc>,7>z)  .  Figure  13  shows  the  power  spectrum  of  two  other  test  signals  sin(co27«)  and 
cos(co2Tn).  Both  sets  of  test  signals  in  these  figures  have  an  SXR  of  20  dB.  The  fre- 
quency/ is  13.45  Hz  and/;  is  23.-45  Hz  and     T  is  — -      second. 

0*4 

Theoretically,  the  power  spectrum  of  these  signals  have  two  impulse  functions  at/ 
and  —  f,  (  /=  1  ,2  ).  Figures  12  and  13  show  experimental  results  agreeing  with  the 
theoretical  results.  Figures  12  and  13  indicate  which  frequency  has  high  power  but  they 
provide  no  phase  information. 

j  The  AR  model  is  sensitive  to  frequency  but  is  not  sensitive  to  the  phase.  For  any 
selection  of  frequency  /  which  satisfies  the  Nyquist  theorem,  the  AR  model  provides 
good  frequency  estimation  provided  the  SXR  is  sufficient  large. 
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Figure   12.       AR  model  performance  test  1  (SNR=  20  dB). 
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Figure   13.        AR  model  performance  test  2  (SNR=  20  dB). 
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B.      DOPPLER  ESTIMATION 

In  section  A  of  chapter  IV,  we  discussed  five  items. 

1.  BPF1  and  BPF2  outputs  have  two  sinusoidal  components  each. 

2.  For  BPF1.  the  dominant  frequency  is  designated  as  wx  and  for  BPF2,  the  dominant 
frequency  is  designated  as  cov  where 

2zlf+  a,) 
-2n(f+  a.-,) 


'*-— Jj 

3.  The  power  spectrum  of  cross  term  of  BPF1  and  BPF2  has  four  components  which 
can  be  detected  in  very  high  SXR  (i.e.  SNR  =  100  dB). 

4.  When  the  SXR  is  low  (i.e.  SXR  =  20  dB).  only  one  dominant  frequency  compo- 
nent is  detected. 

5.  The  dominant  frequency  of  the  power  spectrum  PxyxJJ)  corresponds  to  the  differ- 
ential Doppler  frequency. 

Figure  14  shows  the  power  spectrum  of  the  BPF1  output  when  /  is  13  Hz.  a,  is 
0.45  Hz,  and/  is  64  Hz.  Figure  15  shows  the  power  spectrum  of  the  BPF2  output 
when  /is  13  Hz.  c/.2  is  —0.45  Hz.  and/  is  64  Hz.  As  expected  both  figures  show  two 
spectral  lines  when  the  SXR  is  100  dB.  The  dominant  frequencies  are  located  at/equal 
13.45  Hz  for  BPF1,  and  at /equal  -12.55  Hz  for  BPF2.  Figure  16  shows  the  power 
spectrum  of  the  cross  term  of  BPF1  and  BPF2.  This  figure  shows  four  spectral  lines  at 
an  SXR  of  100  dB.  while  only  one  dominant  frequency  is  detected  at  an  SXR  of  20  dB. 
Dominant  frequencies  are  always  located  between  —1  Hz  and  1  Hz.  Figure  17  shows 
a  subplot  of  Figure  16  for  frequency  between  —1  Hz  and  1  Hz.  This  figure  shows  the 
dominant  frequency /=  .9  Hz  and  the  smaller  spectral  component  at/t  —.9  Hz  for  an 
SXR  of  100  dB.  Only  one  dominant  frequency  can  be  detected  at /equal  0.9  Hz  for  an 
SXR  of  20  dB.  Doppler  difference  from  the  dominant  spectral  component  corresponds 
to  the  true  Doppler  difference  otj  —  a2  which  is  0.9  Hz 

Figure  18  shows  a  comparison  for  SXR  of  20  dB  and  10  dB.  Figure  19  is  subplot 
of  Figure  IS  for  -1  Hz  </<  1  Hz.  From  this  figure,  when  SXR  =  10  dB,  the  dominant 
frequency  is  shifted  a  little  bit  from  the  true  Doppler  difference.  We  have  shown  the  five 
issues  as  we  have  discussed  in  the  previous  chapter.  In  summary,  for  SXRs  greater  than 
or  equal  to  20  dB.  the  power  spectrum  of  the  cross  term  provides  good  Doppler  esti- 
mation. 
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Figure   14.       Power  spectrum  of  the  BPF1  output. 
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Figure   15.       Power  spectrum  of  the  BPF2  output. 
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Figure  16.       Cross  pouer  spectrum  of  BPF1  and  BPF2. 
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Figure  17.       Subplot  of  the  cross  power  spectrum. 
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Figure   18.    .  Comparison  the  Doppler  estimation  at  tuo  different  SNR's. 
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Figure   19.       Subplot  of  the  Doppler  estimation. 
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C.      DIFFERENTIAL  TIME  DELAY  ESTIMATION 

In  Chapter  IV  section  B.l.a  we  discussed  the  problems  associated  with  using  un- 

normalized  quantities  and  suggested  two  possible  approaches  to  solve  the  problem.    In 

this  section  we  will  show  that  for  the  case  of  la,|  <—    Hz  ,  the  power  spectrum  of  a 

6 

nondelayed  signal  is  smaller  than  the  power  spectrum  of  a  signal  delayed  by  one  second 
(  i.e..  see  discussion  following  Eq.  (4.35)  ). 

To  estimate  the  time  delay,  we  suggest  the  two  different  approaches.  The  First  ap- 
proach is  the  normalization  of  the  BPF2  output.  The  second  approach  normalizes  the 
cross  term  Pxyxy{f)  with  the  power  spectrum  of  BPF1  Pxx(j)  and  the  power  spectrum  of 
BPF2  Pyyif)  through  the  equation  defined  by  Eq.  (4.40). 

The  following  figures  are  simulations  for  /=  23  Hz,  a,  =  0.23  Hz,  a2  =  —0.02  Hz, 

T=-—  second,  delay  of  0  and  an  S.XR  of  100  dB.    Figure  20  shows  the  surface  plot  of 
Cm 

Pxyxy(f)  while  Figure  21  shows  the  contour  plot  of  Pxyxy(f).    Figure  22  shows  the  cross 

section  plot  of  Figure  20  with  respect  to  the  DELAY  axis.    This  figure  shows  that  the 

power  spectrum  has  a  spurious  peak  at  delay  of  one  second.    At  the  proper  delay  (i.e., 

0  second  delay)  Pxxxy(f)  has  a  relatively  small  value. 

PxuJf)  is  affected  bv  two  terms,  namelv  pm    and    -— ttt  .  Fisure  23  shows  the 

'  \Axyxv{f)V 

maximum  power  spectrum  term  of  the  transfer  function  and  Figure  24  shows  the  driving 

noise  variance.   At  a  delay  of  ±  1  second  and  of  0  second  the  transfer  function  shows  a 

peak  while  the  variance  of  driving  noise  has  small  value.  In  this  case.  Pxyxy(f),  the  product 

of  the  form  n  ,  — tt—  at  0  second  delav  has  a  smaller  value  than  that  at  any  other 

\Axxxv(f)\2 
delay  location.    Using  the  scheme  shown  in  Figure  10.  we  normalize  the  BPF2  output 

and  present  the  result  in  Figure  25.  We  can  detect  the  proper  time  delay  and  the  proper 
Doppler  difference.  Figures  26  and  27  show  the  corresponding  contour  plot  and  the 
power  spectrum  of  the  normalized  cross  term.  For  any  value  of  time  delay  and  Doppler 
difference,  this  approach  provides  good  information  about  differential  Doppler  and  dif- 
ferential time  delay  provided  the  SXR  is  greater  than  or  equal  to  70  dB. 

Results  of  the  second  approach  (  Figure  11  )  are  demonstrated  in  Figures  28  and 
29.  This  approach  also  gives  good  information  about  differential  time  delay  and  differ- 
ential Doppler.  This  approach  provides  good  time  delay  and  Doppler  estimation  for 
SXR  greater  than  or  equal  to  20  dB. 

If  we  use  the  contour  shape,  we  can  also  estimate  the  differential  time  delay  and 
differential  Doppler  for  SXR  of  less  than  20  dB.  Differential  time  delay  and  differential 
Doppler  estimation  can  be  extended  down  to  SXRs  of  about  0  dB  when  using  the  con- 
tour of  the  coherence  surface.   Figure  30  is  a  contour  plot  of  using  this  second  approach 
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(  i.e..  coherence  approach  )  for  a,  =  0.23  Hz,  a2  =  —  0.02  Hz,  Delay  =  0,  and  an  SXR  of 
0  dB.  The  contour  shape  of  the  second  approach  has  the  shape  of  the  letter  X.  When 
the  SXR  is  greater  or  equal  to  0  dB.  the  proper  time  delay  and  Doppler  difference  seems 
to  be  located  at  the  cross  over  point  of  the  X.  When  the  SXR  is  less  than  0  dB,  the 
contour  shape  does  not  shows  an  X  clearly.  In  the  coherence  approach,  Figure  31 
provides  the  contour  plot,  with  channel- 1  containing  only  noise  and  channel-2  contain- 
ing signal  plus  noise  (  SXR  =  0  dB  ).  Figure  32  shows  the  contour  plot,  when  channel-2 
contains  only  noise  and  channel- 1  contains  signal  plus  noise  (  SXR  =  0  dB  ). 
Figure  33  shows  the  contour  plot,  when  both  channels  contain  noise  only.  The  above 
three  figures  do  not  show  an  X  clearly,  hence  do  not  allow  signal  detection  for  esti- 
mation of  any  parameter.  Using  this  information,  we  can  distinguish  the  signal  combi- 
nations from  the  noise  combinations. 

When  we  estimate  the  time  delay  and  Doppler  difference  using  the  contour  plot,  the 
modified  cross  power  spectral  approach  has  some  problems  as  discussed  in 
Appendix  D. 
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Figure  20.       Surface  plot  of  the  power  spectr 
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Figure  21.       Contour  plot  of  the  power  spectrum. 
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Figure  22.        Maximum  power  spectrum. 
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Figure  23.       Maximum  power  spectrum  of  the  transfer  function. 
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Figure  2-4.       Variance  of  the  driving  noise. 
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Figure  25.       Maximum  modified  cross  power  spectrum  of  the  transfer  function. 
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Figure  26.        Power  spectrum  of  the  transfer  function  (surface). 
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Figure  27. 


Power  spectrum  of  the  transfer  function  (contour). 
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Figure  28.       Surface  plot  of  the  coherence. 
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Figure  29.        Contour  plot  of  the  coherence. 
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Figure  30.        Estimation  using  the  contour  shape  (low  SNR). 
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Figure  31.       Contour  plot  (noise  in  channel- 1  only). 
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Figure  32.       Contour  plot  (noise  in  channel-2  only). 
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Figure  33.       Contour  plot  (noise  only). 
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VI.      CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  the  differential  time  delay  and  differential  Doppler  are  estimated  using 
two  different  approaches.  The  first  approach  uses  a  modified  cross  power  spectrum  al- 
gorithm while  the  second  approach  uses  a  coherence  algorithm.  In  both  cases 
autoregressive  modeling  of  the  required  cross  spectral  component  is  used.  The  differen- 
tial time  delay  and  Doppler  difference  can  be  estimated  using  a  threshold  at  high  S.XRs 
and  a  contour  plot  at  low  SXRs. 

At  high  SXRs  (  i.e.,  SXR  >  70  dB  ).  the  first  approach  locates  the  dominant  peak 
at  the  proper  time  delay  and  the  proper  Doppler.  The  second  approach  utilizes  two 
additional  AR  models  and  two  additional  FFTs  to  obtain  the  autopower  spectrum  of 
channel- 1  and  channel-2.  At  moderate  SXRs  (  i.e.,  SXR  >  20  dB  ).  the  second  approach 
has  the  highest  coherence  peak  at  the  proper  time  delay  and  Doppler  frequency.  When 
the  SXR  >  0  dB.  the  contour  plot  of  the  coherence  has  the  shape  of  the  letter  X.  The 
differential  time  delay  and  Doppler  difference  can  be  located  at  the  cross  over  point  of 
the  X.  In  high  SXRs  we  can  use  the  highest  peak  to  detect  the  time  delay  and  Doppler 
frequency  using  either  approach.  When  we  use  the  contour  plot,  the  information  about 
time  delay  and  Doppler  depends  on  the  location  and  the  form  of  the  cross  over  point 
of  the  X. 

We  can  estimate  the  differential  time  delay  and  differential  Doppler  using  the  AR 
modeling  of  coherence.  But  we  can  not  get  the  numerically  correct  value  of  the  coher- 
ence coefficient.   The  following  three  points  are  recommended  for  future  study. 

1.  To  get  a  cross  power  spectrum,  we  assumed  that  |  Pxyxy{f)  |  ~  |  Pxy{f)  | 2  .  This  was 
required  since  we  can  not  get  a  cross  power  spectrum  using  an  AR  model.  Im- 
provement of  the  cross  power  spectrum  estimates  using  an  AR  model  can  lead  to 
improvement  of  the  algorithm. 

2.  Our  first  approach  is  not  normalized  by  the  auto  power  spectra.  If  we  find  the 
corresponding  auto  power  spectra,  we  can  also  obtain  a  coherence  coefficient  esti- 
mate. 

3.  To  get  a  better  information  at  a  low  SXR,  we  need  to  study  the  characteristics  of 
the  contour  plots. 


54 


APPENDIX  A.      PHASE  DERIVATION 

Define  the  phase  at  instant  k. 


2n(f+oix)k 
A=         fs  (4-4) 

2n<f+  a2)k 

A=    \  (4-5) 


Let  x{n)  denote  the  sampled  analog  signal  x{t)  \  !=nT 
then 


x{n)  =  cos{2n{f+  a,)  —  } 


Similiarilv 


=  cos{2tt(/'+  a,)  —  +  >0j 


*(«  +  l)  =  cos{27r(/'+a1)-!17:J-} 

=  cos{2^+^)7-+       '    ■     ]     } 
=  cos(27r(/+a1)-^-+^1} 


*(«  +  2)  =  cos{2n(f+  aj)  -^7^"  } 


fa.!' 


=  cos{2n(f+  a,)  f-  +       V  .    "     }  («.3] 

=  cos{27rOr+a1)7-  +  ^2} 


x(«  +  m)  =  cos{2n(f+  a,)  — - —  } 

Js 

n        2n(f+  a,)m 
=  cos{2*(f +  a,)  f  +  '        }  (4.7) 

=  cos{2n(f+<xl)f  +  x<t>m} 

Js 
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y{n)  =  cos{2n{f+  a2)  -j- ) 

Js 

=  cos{2t(/'+  a2)  -j-  +  yd> 

Js 


(a.5) 


y{n  +  1)  =  cos{2n(j'+  a2)  ~^-  } 

Js 

n        2-{f+  a2) 
=  cos{2tt(/"+  a2)  —  + — —  } 

=  cos{2n{f+cc2)-j-  +  y<{)2} 
Js 


(a.6) 


j-(/7  +  2)  =  cos{27rf/+a2)^^-} 

„        2n{f+  a2)2 

=  cos{2;:(7+  a2)  -5-  +  — v        2'     } 

Js  Js 

-cos{2ic(/,+  a2)^-  +  ^2} 


(a.7) 


>•(«  +  m)  =  cos{2tt(/*+  a2)  -^"L } 

,,        2n(f+  a2)w 
=  cos{27r(/'+a2)^-+       U 

=  cos{2;rl/+c<2)-JL  +  v(/»w} 


(4.8) 
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APPENDIX  B.      OUTPUT  OF  BPF2 

Define  two  functions  depending  on  k. 


A--) -A: 


kAy=    y  cos(2;r/;-f  y'27r/A-  (4.25) 


«=o 


N~\-k 

kBy  =    £  M^fif)^1^  (4.26) 


where      f2=f+  »;• 

A.      FULL  INFORMATION  CASE  (  N  POINTS  ) 

A'  points  are  data  points  at  BPF2  input. 


.Y-l 

Y*JJ)  =  Y  cos(2^2  f  +  ^J^'f 

fV-1 


En  n  i2-f— 

{  cos(2tt/2  —  )  cosy<}>m  -  sm{2nf2  — )  sm/J.f'  '    <v 

A'-l  A'-l 


=  cos^^^cosl^^y-)^'2"  v  -  sin>,0m2^sin(2^/2y-)c  ;2r:  a 
=  0AV  cos  V<j)m-0BV  sin  ..</>, 


0^>'  lw^m  —  O^y  511i>^m 


/v<*m  +  ^-A<?m  ^yd>„  _  e-jy4>r, 


-oA  2  -0By 

0Ay+j0By  0Av-j0Bv 

2  2 
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B.      PARTIAL  INFORMATION  CASE  (  N-l  POINTS  ) 

iV  —  1  points  are  data  points  at  BPF2  input. 


.V-l-l 


Y*m(f)=     £  cos{2nf2jr  +  y<i>m)^fnN 


73=0 

A"- 1-1 


=    )     {  cos(2rzf2  j- )  cos y4>m  -  sm{2nf2  -j  )  siny/vie72"7  «v 

A -1-1  A- 1-1 

=  cos y<f>m   )^  cos(2^/;y-)^2r:~-sin>,(/)OT  2^  sin(27:/2  -j  )e~j2*  ~ 
=  xA,cosy(pm-xB,sm,4>m 

=  1  A  i  -  1  By  ~ 


(b.2) 


2/ 


AV+j\By  lAy-j.By 

: _  g/y*  m   _|_  : •_  e     Jy^r, 


C.      PARTIAL  INFORMATION  CASE  (  N-K  POINTS  ) 

A'—  k  points  are  data  points  at  BPF2  input. 


A-i-l: 

}-;(/)=  T  cos(27if2f+y<f>„y2*fJk 

.y-i-fc 
=    )     {  cos(2-y;  —  )  cos  y<f>m  -  sin(2re/2  — )  sin y4>my  *J  x 


/       l  ^s^/c/2     ,.    /^5}vm  —  »uivi'y2    <~ 

n=0 

iV-l-fc  A'-l-A- 

V/7         /2-  —  V^  A7         —ills  — 

_^  cos(2w/2-T-y  "  .v  -siny/>m  2^  sin(27r/2  —  )e  ^"  .v      ^3] 

n=0  /j  »=0  /j 


-  kAy  ]  ~  *    >'  2/ 


kAy+jkBy  kAv-jkBv 

— — e>  m  +  — — ; e  Jy  " 
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APPENDIX  C.      AMPLITUDE  COMPARISON 

Let  C,  and  C2  be  the  complex  amplitude. 

1  -  ^ 


Q  = 


/2™47 


1  -  e12- 


c2  = 


—  e 


x_e-Jl*Mr4lf 


(c.l) 


(c.2) 


where  |  a  I  <  0.5,  and  N  >  2f+  1  >  2(/+  a).  The  numerators  of  I  C,  I  and  |  C2 1  are  same, 
because  the  magnitude  of  the  complex  conjugate  is  the  same.  If  we  show  that 
I  1  -  eJ2n'  a7  |  <  |  1  -  e-j2n{2'")'s\  ,  then  we  can  say  j  C, |  >  \C2\.  As  shown  in  Figure  3-4 
|  i  _  e*«*-k\  =  1.-1/ 1  and  1 1  -  g-;W>i|  =  \Bl\  . 


1 

Im 

\a 
^^f9i    \i               Re 

1                o 

xJe2  j 

J/B               91-2TT0i/N 
— — ^                   62=  -2TT(2f+(*)/N 

Figure  34.       Magnitudes  of  two  complex  number. 


2n(2f+a)         2n(N-2f-u)        2n(\  -  a) 
The  smallest  anele  LIOB  is — or — =  ■ 


For  positive/ 


L10B  = 


X 


2n(2f+  a) 

N 


N 


> 


Inty. 
N 


=  LlOA 


(c.3) 
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Since  |  1  -  y.  |  >  |  a 


LIOB  = 


2n(\  -a) 
A' 


> 


A" 


=  LIOA 


From  Eq.  (c.3)  and  Eq.  (c.4)  \AI\  <\BI\  ,  therefore 

lc.|>lc,| 


(c.4) 


(c.5) 
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APPENDIX  D.     CONTOUR  PLOTS  OF  MODIFIED  CROSS  POWER 

SPECTRUM 

For  SXRs  greater  than  or  equal  to  10  dB  the  contour  plot  of  the  modified  cross 
power  spectrum  algorithm  has  the  shape  of  a  boomerang.  The  shape  of  the  boomerang 
depends  on  a,.  Figure  35  and  Figure  36  show  the  contour  plot  of  the  modified  cross 
power  spectrum.  The  two  figures  show  the  boomerang  shape  and  their  different  di- 
rections. For  positive  a2  (  case  1  )  and  negative  a2  (  case  2  )  the  opening  of  the 
boomerang  is  toward  the  left  and  right,  respectively.  In  the  both  cases,  we  can  estimate 
the  time  delay  and  Doppler  frequency  using  these  contour  plots.  The  proper  time  delay 
and  Doppler  difference  can  be  detected  by  locating  the  point  of  symmetry  of  the 
boomerang.  The  Figure  37  shows  a  contour  plot  when  a2  is  —0.02  Hz.  Since  a2  is  very 
small,  the  contour  plot  looks  like  a  straight  line.  In  this  case,  the  differential  time  delay 
and  Doppler  difference  can  still  be  estimated  by  locating  the  point  of  symmetry  of  the 
boomerang  (  i.e.,  center  point  ).  Figure  3S  shows  the  contour  plot,  when  channel-1 
contains  noise  only  and  channel-2  contains  signal  plus  noise  (  i.e..  SXR  =  10  dB  ).  Fig- 
ure 39  shows  the  contour  plot,  when  channel-2  contains  noise  only  and  channel-1  con- 
tains signal  and  noise  {  i.e..  SXR  =  10  dB  ).  Figure  40  shows  the  contour  plot,  when 
both  channels  contain  noise  only.  Using  the  modified  cross  power  spectrum  technique, 
the  three  types  of  noise  only  set  ups  appear  as  similiar  plots  (  see  Figures  38.  39  and 
40  ).  Figure  38  seems  to  indicate  the  presence  of  signals  in  both  channels,  while  Figures 
39  and  40  appear  more  like  noise.  Hence,  contour  plots  of  the  modified  cross  power 
spectrum  do  not  allow  clear  identification.  In  this  case,  it  will  be  difficult  to  distinguish 
between  the  noise  only  situation  and  a  small  a2  situation. 

To  use  the  modified  cross  power  spectrum  algorithm  at  low  SXRs  (  i.e.. 
SXR  <  70  ).  we  need  to  investigate  the  contour  plot  some  more. 
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Figure  35.       Contour  plot  of  the  modified  cross  power  spectrum  (case  1). 
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Fi»ure  36. 


Contour  plot  of  the  modified  cross  power  spectrum  (case  2). 
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Figure  37.       Contour  plot  of  the  modified  cross  power  spectrum  (case  3). 
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Figure  38.        Contour  plot  of  the  modified  cross  power  spectrum  (noise  in  channel- 1 
only). 
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Figure  39.        Contour  plot  of  the  modified  cross  power  spectrum  (noise  in  channel-2 
only). 
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Figure  40.       Contour  plot  of  the  modified  cross  power  spectrum  (noise  only). 
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APPENDIX  E.      MODIFIED  CROSS  POWER  SPECTRUM  PROGRAM 

J*  -'.  ^'.  .'»  »'-  .',  .'-,  »•-  »i-  y-  *i.  »'*  »'-  w*-  -'-  -»-  ..».  ..'-  Ji.  *)-  >.'-  ij'-  ..'-  *'-  *'*  •'*  »«-  **-  »'-  y-  .J-  .J*  y,  *«.  y^  ^  ^  »*-  y-  .J-  J-  .'.  J-  .J-  J,  »'.  »'*  J*  J*  »'*  »»*  JL  j*  „'.  J«  -'-  »J-  J-  «J-  J-  JL  -J-  -'..  -!-  -J-  -*-  J-  -'.. 

c 

c      main  program  of  approach  one.  * 

c  this  program  share  subroutine  with  * 

c  the  main  program  approach  two  * 

c 

COMPLEX  X(-128: 255),Y(-128: 255), A(0:  2500), B(0:  2500) 

COMPLEX  BPFK-128:  191)  ,BPF2( -128:  191),BPF(0:  100) 

INTEGER  SNR,BINDEL, DELAY 

DATA  X,Y,A,B/5770*(0. ,0. ) / ,BPF1 ,BPF2 ,BPF/741*( 0. ,0. )/ 

0PEN(UNIT=7,FILE=' FILENAME  FILETYPE  A') 

PI=AC0S(-1.  ) 

c - - * 

c 

c     input  part  * 

c  * 

c     fs      -  sampling  frequency  * 

c     f       -  carrier  frequency  * 

c     alphal   -  doppler  effect  at  sensorl 

c     alpha2   -  doppler  effect  at  sensor2  * 

c     delay    -  signal  receiving  time  difference  * 

c  between  two  sensors  * 

c     SNR     -  signal  to  noise  ratio  * 

c     iseed    -  noise  generator  seed(odd  number)  * 

en       -  number  of  data  during  one  second  - 

c 

c * 

FS=64. 
F=23. 

ALPHA1=. 23 
ALPHA2=-.  02 
DELAY=0 
SNR=100 
N=64 

T=2.  *PI*F/FS 
ISEED=13 
c * 

c 

c        input  data  sampling  at  two  sensors  * 

c  * 

c * 

CALL  SIGNAL(X,F+ALPHA1,FS,0, SNR, ISEED, 0,0) 
CALL  SIGNAL(Y,F+ALPHA2,FS, DELAY, SNR, ISEED, 0,0) 
DO  100  K=-128,191 
DO  50  J=K,K+63 


6S 


c ,r 

C  " 

c  linear  transformation  using  band  pass  filter 

c  * 

c * 

OMEGA=T*FLOAT(J-K) 

BPF1(K)=BPF1(K)+X(J)*CMPLX(C0S(0MEGA),-SIN( OMEGA)) 

50  BPF2(K)=BPF2(K)+Y(J)*CMPLX(C0S(0MEGA),SIN(0MEGA)) 

c * 

c  * 

c  normalization  of  BPF2  * 

c 

c * 

100       BPF2(K)=BPF2(K)/CABS(BPF2(K)) 

DO  300  BINDEL=-128,127 

DO  200  K=0,N 

c - * 

c 

c  input  of  AR  model  * 

c 

c -,v 

200  BPF(K)=BPF1(K)*BPF2(K+BINDEL) 

CALL  CLEAR(A,B,2500) 
CALL  BURGAR(BPF,N+1,A,IP,VAR) 
CALL  CHANGEFFT(B,A,11,MAX) 
DO  250  I=-32,-l 
250      V,"RITE(7,997)BINDEL,FLOAT(I)/32.  ,CABS(1.  /B( 2048+1 ) )*VART 
*,CABS(1. /B(2048+I)) 
DO  300  1=0,32 
300      WRITE(7,997)BINDEL,FLOAT(I)/32.  ,CABS( 1. /B( I) )*VART 
*,CABS(1.  /B(I)) 
CL0SE(7) 
STOP 

997    F0RMAT(1X,I4,2X,E13.  6,2X,E13.  6,2X,E13.  6) 

END 


69 


APPENDIX  F.      COHERENCE  PROGRAM 

»t-  JL  -T-  jl  jl  JL  JL  J-  JL  jl  y-  jl  jl  jl  jl  -i-  jl  jl  jl  y-  a  jl  JL  jl  jl  jl  jl  -f-  J.  J/f  J-  jl  JL  -»-  «'-  JL  J*  jl  jl  »'f  y-  J-  J*  j*  jl  j-  jl  jl  j-  jl  jl  jl  jl  J-  jl  jf  jl  jl  jl  jl  j-  jl  jl  jl  jl  jl  jl  jl.  jl 

c  * 

c     main  program  of  approach  two  * 

c  * 

JL  J-  -'-  JL  JL  -•-  JL  JL  Jl-  *'-  -■-  ^-  J-  JL  JL  JL  JL  JL  JL  J-  JL  JL  JL  JL  JL  JL  JL  JL  J/-  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  J-  JL  JL  -L  J-  JL  JL  J-  JL  JL  JL  JL  JL  J-  JL  JL  JL  JL  JL  JL  J*.  JL  JL  JL  JL  JL 

COMPLEX  X(-128:  255),Y(-128: 255), A(0: 2500), B(0: 2500) 

COMPLEX  BPFK-128:  191)  ,BPF2( -128:  191),BPF(0:  100) 

COMPLEX  AUTOXC0: 100) , AUT0Y( 0: 100),SUMXY 

INTEGER  SNR,BINDEL, DELAY 

DATA  X,Y,A,B/5770*(0. ,0. ) / ,BPF1 ,BPF2 ,BPF/741*( 0. ,0. )/ 

DATA  AUTOX,AUTOY/202*(0. ,0. )/ 

0PEN(UNIT=7,FILE=' FILENAME  FILETYPE  A') 

PI=AC0S(-1.  ) 

c * 

c  * 

c     input  part  * 

c  * 

c     fs      -  sampling  frequency 

c     f       -  carrier  frequency  * 

c     alphal   -  doppler  effect  at  sensorl  * 

c     alpha2   -  doppler  effect  at  sensor2  * 

c     delay   -  signal  receiving  time  difference  * 

c  between  two  sensors 

c     SN'R     -  signal  to  noise  ratio 

c     iseed   -  noise  generator  seed(odd  number)  * 

en       -  number  of  data  during  one  second  * 

c  * 

c ft 

ALPHA1=. 23 

ALPHA2=-.  02 

SNR=20 

F=23. 

N=64 

FS=64. 

T=2. *PI*F/FS 

ISEED=23 

DELAY=0 

c * 

c  * 

c     input  data  sampling  at  two  sensors  * 

c  * 

c * 


CALL  SIGNAL(X,F+ALPHA1,FS,0,SNR, ISEED, 0,0) 
CALL  S IGNAL( Y , F+ALPHA2 , FS , DELAY , SNR , ISEED ,1,0) 
DO  100  K=-128,191 
DO  50  J=K,K+63 


c  * 

c     linear  transformation  using  band  pass  filter  * 

c 


c- 
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OMEGA=T*FLOAT(J-K) 

BPF1(K)=BPF1(K)+X(J)*CMPLX(C0S(0MEGA),-SIN(0MEGA)) 
50  BPF2(K)=BPF2(K)+Y(J)*CMPLX(COS(OMEGA),SIN(OMEGA)) 

100   CONTINUE 
SUMX=0. 
DO  150  K=0,N 

c * 

c 

c     find  AR  order  at  sensorl  * 


150      AUT0X(K)=BPF1(K) 

CALL  CLEAR(A,B,2500) 

CALL  BURGAR(AUT0X,N+1,A,IP,VARX) 

CALL  CHANGEFFT(B,A,11,MAX) 
XMAX=CABS(1. /B(MAX)) 
DO  300  BINDEL=-128,127 

DO  200  K=0,N 

c * 

c  * 

c     input  of  AR  model  * 

c 

c -v 

AUT0Y(K)=BPF2(K+BINDEL) 

200  BPF(K)=BPF1(K)*(BPF2(K+BINDEL)) 

c -,v 

c  * 

c       find  auto  power  spectrum  of  sensor2  at  delay  k  * 

c 

c * 

CALL  CLEAR(A,B,2500) 

CALL  BURG(AUT0Y,N+1,A,IP,VARY) 

CALL  CHANGEFFT(B,A,U,MAX) 

YMAX=1. /CABS(B(NAX)) 

c Vr 

C  " 

c       find  auto  power  spectrum  of  cross  term  * 

c  * 

c -v 

CALL  CLEAR(A,B,2500) 
CALL  BURG(BPF,N+1,A,IP,VARXY) 
CALL  CHANGEFFT(B,A,11,MAX) 
XSMAX=(XMAX*YMAX)**2*VARX*VARY/VARXY*fs 
DO  250  I=-32,-l 
250   WRITE(7,997)BINDEL,FLOAT(I)/32. ,CABS(1. /B(2048+I) )/SQRT(XSMAX)*8. 
'V,CABS(1.  /B(2048+I)) 
DO  300  1=0,32 
300      WRITE(7,997)BINDEL,FLOAT(I)/32.  ,CABS( 1. /B( I) )/SQRT(XSMAX)*8. 
*,CABS(1.  /B(I)) 
CL0SE(7) 
STOP 

997   F0RMAT(1X,I4,2X,E13.  6,2X,E13.  6,2X,E13.  6) 

END 
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SIGNAL  GENERATER 
INPUT 


OUTPUT 


N      -  NUMBER  OF  DATA  POINTS 
F1,F2   -  FREQUENCY  OF  1ST, 2ND  SIGNAL 
AMP 1,AMP2 -FREQUENCY  OF  1ST, 2ND  SIGNAL 


A(N) 


SIGANL 


* 

i'- 

Vr 
* 

* 


SUBROUTINE  SIGNAL(B ,F,FS ,D, SNR, ISEED, ID, NOISE) 

COMPLEX  B(-128: 255) 

REAL  A(0:  1) 

INTEGER  SNR,D 

PI=ACOS(-l.  ) 

DATA  A/0.  ,0.  / 

SIGMA=1. 

A(0)=SQRT(2.*10.**(SNR/10)) 

IF(NOISE.EQ.  1)  A(0)=0. 
DO  100  I=-12S,255 

M=(I-D)/128 

IF(I.  LT.  D.  OR.  M.NE.  0)  M=l 

T=AMOD(F*FLOAT( I-D) ,FS) 

CALL  GAUSS( ISEED, SIGMA, 0.  , RANDOM) 

IF(ID.EQ.  0)  THEN 

X=COS(2.*PI*T/FS)*A(M)+RANDOM 

ELSE 

X=SIN(2.*PI*T/FS)*A(M)+RAND0M 

END  IF 
100      B(I)=CMPLX(X,0.  ) 
RETURN 
END 


BURG  ALGORITHM 
INPUT 


OUTPUT 


N      -  NUMBER  OF  DATA  POINTS 
X      -  INPUT  SIGNAL 

IP     -  ORDER  OF  AR 

A(0:IP)-  AR  COEFFICIENTS 

VAR    -  DRIVING  NOISE  VARIATION 


■>.- 
* 

* 
* 


SUBROUTINE  BURGAR(X,N, A, IP, VAR) 

COMPLEX  X(N) ,A(0:  N) ,EFK( 500) ,EBK(500) 

COMPLEX  EFK1(500),EBK1(500),AA(20,20),SUMN,SUMD 

REAL  RHO(0: 80),FPE(0: 80) 

INTEGER  START 

RHO(0)=0. 

FPE(0)  =  1.E64 

A(0)=CMPLX(1.  ,0.  ) 

IP=1 
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START=1 
DO  10  1=1, N 
10       RHO(0)=RHO(0)+CABS(X(I))**2/FLOAT(N) 
DO  20  1=2, N 

EFK1(I)=X(I) 
20       EBK1(I-1)=X(I-1) 
LOOP 
K=IP 

SUMN=CMPLX(0. ,0. ) 
SUMD=CMPLX(0. ,0. ) 
DO  30  I=K+1,N 

SUMN=SUMN+EFK1(I)*C0NJG(EBK1(I-1)) 
30  SUMD=SUMD+CABS(EFK1(I)**2)+CABS(EBK1(I-1)**2) 

CO  SUMD=SUMD+CABS(EFK1(I))**2+CABS(EBK1(I-1))**2 

AA(K,K)=-2. *SUMN/SUMD 
TEMP=1.  -CABS(AA(K,K))**2 

IFCTEMP.  LE.  0.  )  TEMP=1.  E-10 
RH0(K)=TEMP*RH0(K-1) 
IF(K. GT.  1)  THEN 
DO  40  J=1,K-1 
40  AA(J,K)=AA(J,K-1)+AA(K,K)*C0NJG(AA(K-J,K-1)) 

ENDIF 
DO  60  I=K+2,N 

EFK(I)=EFK1(I)+AA(K,K)*EBK1(I-1) 
60  EBK(I-l)=EBKl(I-2)+C0NJG(AA(K,K))*EFKl(I-l) 

DO  70  I=K+2,N 

EFK1(I)=EFK(I) 
70       EBK1(I-1)=EBK(I-1) 

IF(N-K.  EQ.  1)  THEN 

FPE(K)=FPE(K-1)+1. 
ELSE 

FPE(K)=RH0(K)*FL0AT(N+1+K)/FL0AT(N-1-K) 
ENDIF 
IP=IP+1 
UNTIL(  FPE(K).  GT.  FPE(K-l) .  AND.  K.  GT.  START) 
IP=K-1 

DO  100  1=1, IP 
100      A(I)=AA(I,IP) 
VAR=RHO(IP) 
RETURN 
END 


C * 

c  * 

C     BURG  ALGORITHM 

C  INPUT 

C  N        NUMBER  OF  DATA  POINTS 

C  X        INPUT  SIGNAL  * 

C  IP       ORDER  OF  AR 

C  OUTPUT 

C  A(0:IP)-  AR  COEFFICIENTS 

C  VAR    -  DRIVING  NOISE  VARIATION  * 

C  * 


C- 


73 


SUBROUTINE  BURG(X ,N, A, IP , VAR) 
COMPLEX  X(N) ,A(0:  N) ,EFK(500) ,EBK(500) 
COMPLEX  EFK1(500) ,EBK1(500) ,AA( 20,20) ,SUMN,SUMD 
REAL  RHO(0: 80),FPE(0: 80) 
INTEGER  START 
RHO(0)=0. 
A(0)=CMPLX(1.  ,0.  ) 
DO  10  1=1, N 
10       RHO(0)=RHO(0)+CABS(X(I))**2/FLOAT(N) 
DO  20  1=2, N 

EFK1(I)=X(I) 
20       EBK1(I-1)=X(I-1) 
DO  1000  K=1,IP 

SUMN=CMPLX(0. ,0. ) 
SUMD=CMPLX(0.  ,0.  ) 
DO  30  I=K+1,N 

SUMN=SUMN+EFK1( I )*CONJG( EBK1(  I - 1 ) ) 
30  SUMD=SUMD+CABS(EFK1( I) )**2+CABS(EBKl( 1-1) )**2 

AA(K,K)  =  -2.  '-SUMN/SUMD 
TEMP=1.  -CABS(AA(K,K))**2 

IF(TEMP.  LE.  0.  )  TEMP=1.E-10 
RHO(  K)=TEMP'»-RHO(  K- 1 ) 
IF(K.  GT.  1)  THEN 
DO  40  J=1,K-1 
40  AA(J,K)=AA(J,K-l)+AA(K,K)*CONJG(AA(K-J,K-l)) 

END  IF 
DO  60  I=K+2,N 

EFK(I)=EFK1(I)+AA(K,K)*EBK1(I-1) 
60  EBK(I-l)=EBKl(I-2)+CONJG(AA(K,K))*EFKl(I-l) 

DO  70  I=K+2,N 

EFK1(I)=EFK(I) 
70       EBK1(I-1)=EBK(I-1) 
1000   CONTINUE 

DO  100  1=1, IP 
100      A(I)=AA(I,IP) 
VAR=RHO(IP) 
RETURN 
END 


C * 

c  * 

C     CHANGE  FFT  * 

C  INPUT 

C  B      -  AR  COEFFICIENTS  * 

C  A        TEMPORARY  USING  ARRAY 

C  ISIZE   -  ALOG(FFT  DATA  POINTS)/ALOG2 

C  OUTPUT  * 

C  A      -  POWER  SPECTRUM 

C  * 

C * 
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SUBROUTINE  CHANGEFFT( A, B , ISIZE ,K) 
COMPLEX  A(0:  2**ISIZE-1)  ,B(0:  2**ISIZE-1) 

INTEGER  REVRSE 
PI=ACOS(-l.  ) 
N=2**ISIZE-1 
K=0 
FS=FLOAT(N+l) 

CALL  FFT(N,ISIZE,A,B,FS) 

CALL  FINDMAX(A,N,K) 
RETURN 
END 

C * 

c  * 

C     FINDMAX 

C  INPUT 

C  N        NUMBER  OF  SPECTRUM  POINTS(N+l)  * 

C  A        POWER  SPECTRUM  * 

C  OUTPUT 

C  MAX    -  ARRAY  INDEX  OF  MAXIMUM  POWER  SPECTRUM 

C 

C * 

SUBROUTINE  FINDMAX ( A,N,MAX) 

COMPLEX  A(0:N) 

MAX=0 

AMAX=1.E66 

DO  100  1=0, N 

IF(CABS(A(I)).LT.  AMAX)  THEN 
MAX=I 

AMAX=CABS(A(I)) 
ENDIF 
100   CONTINUE 
RETURN 
END 

C * 

c 

C  REVERSE 

C  BIT  REVERSE  ORDER  INDEX  CHANGING  SUBROUTINE             * 

C  INPUT                                               * 

C  N      -  ARRAY  INDEX 

C  ISIZE   -  ALOG(FFT  DATA  POINTS )/AL0G2 

C  OUTPUT                                              * 

C  REVERSE-  BIT  REVERSE  ORDER  INDEX                * 

C  * 

C * 

INTEGER  FUNCTION  REVRSE(N, ISIZE) 
INTEGER  A(20) 
DO  200  1=1, ISIZE 
A(I)=MOD(N,2) 

N=N/2 
200   CONTINUE 

REVRSE=0 
DO  300  1=1, ISIZE 
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REVRSE=REVRSE+A(  I  )*2< 
300   CONTINUE 
RETURN 
END 


•(ISIZE-I) 


DFT 


COOLEY  -  TUKEY  ALGORITHM 
INPUT 

N 

ISIZE 
A 
B 


OUTPUT 


NUMBER  OF  DATA  POINTS 

ALOG(N)/ALOG2 

BIT  REVERSE  ORDER  INDEXED  TIME  DOMAIN 

TEMPORARY  ARRAY 


FREQUENCY  DOMAIN 


* 
* 
* 

* 

Vc 
Vc 
* 


200 


300 


400 
500 


SUBROUTINE  DFT(N, ISIZE , A,B) 
COMPLEX  A(0:N)  ,B(0:N),W 
PI=ACOS(-l.  ) 
DO  500  11=1, ISIZE 
ISAGE1=2**(I1-1) 
ISTAGE=2**I1 
DO  200  1=0, N 

ITEST=MOD(I,ISTAGE) 
IF(ITEST. GT.  ISAGE1)  THEN 
K=ITEST-ISAGE1 

THEATA=2. *PI*FLOAT(K)/FLOAT( ISTAGE) 

T1=SIN(THEATA) 

T2=COS(THEATA) 

IF(ABS(T1).LT.  1E-5)  T1=0 
IF(ABS(T2). LT. 1E-5)  T2=0 
W=CMPLX(T2,-T1) 
A(I)=A(I)*W 
END  IF 
CONTINUE 
DO  300  1=0, N 

ITEST=MOD( I, ISTAGE) 
IF(ITEST. GE. ISAGE1)  THEN 

B(I)=-A(I)+A(I-ISAGE1) 
ELSE 

B(I)=A(I)+A(I+ISAGE1) 
END  IF 
CONTINUE 
DO  400  1=0, N 

A(I)=B(I) 
CONTINUE 
CONTINUE 
RETURN 
END 
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FFT 


1. INDEX  CHANGE 

2. USE  COOLEY  TUKEY  ALGORITHM 
INPUT 

N 

ISIZE 

B 

FS 


NUMBER  OF  DATA  POINTS 
ALOG(N)/ALOG2 
INPUT(TIME  DOMAIN) 
SAMPLING  FREQUENCY 


OUTPUT 


FREQUENCY  DOMAIN 


* 
* 
* 

* 
it 

it 
it 
it 


SUBROUTINE  FFT(N, ISIZE , A,B ,FS) 

COMPLEX  A(0:N),B(0:N) 

INTEGER  REVRSE 

DO  100  I=0,N 

K=I 
100      A(REVRSE(K,ISIZE))=B(I) 
CALL  DFT(N,ISIZE,A,B) 

DO  200  I=0SN 
200      A(I)=A(I)/FS 

RETURN 

END 


CLEAR 


INPUT 


OUTPUT 


N      -  NUMBER  OF  DATA  POINTS 
A,B    -  INPUT  arrays 

A,B    -  reinitiallized  arrays 


•it 
it 

V- 
i- 
-V 


100 


SUBROUTINE  CLEAR(A,B,N) 
COMPLEX  A(0:N),B(0:N) 
DO  100  1=0, N 

A(I)=B(I)=CMPLX(0. ,0. ) 
RETURN 
END 
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